Towards adiabatic waveforms for inspiral into Kerr black holes: 
I. A new model of the source for the time-domain perturbation equation 
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We revisit the problem of the emission of gravitational waves from a test mass orbiting and thus 
perturbing a Kerr black hole. The source term of the Teukolsky perturbation equation contains a 
Dirac delta function which represents a point particle. We present a technique to effectively model 
the delta function and its derivatives using as few as four points on a numerical grid. The source 
term is then incorporated into a code that evolves the Teukolsky equation in the time domain as 
a (2+1) dimensional PDE. The waveforms and energy fluxes are extracted far from the black hole. 
Our comparisons with earlier work show an order of magnitude gain in performance (speed) and 
numerical errors less than 1% for a large fraction of parameter space. As a first application of this 
code, we analyze the effect of finite extraction radius on the energy fluxes. This paper is the first 
in a series whose goal is to develop adiabatic waveforms describing the inspiral of a small compact 
body into a massive Kerr black hole. 

PACS numbers: 04.25.Nx, 04.30.Db, 04.30.-w 



I. INTRODUCTION 
A. Background 

The extreme mass ratio limit of binary systems — bi- 
naries with one mass far smaller than the other — has 
been a special focus of research in gravitation in recent 
years. This is in part because this problem is, at least 
formally, particularly clean and beautiful: the mass ra- 
tio allows us to treat the binary as an exact black hole 
solution plus a perturbation due to the secondary mass. 
Perturbative techniques can be used to analyze the sys- 
tem, making it (in principle at least) much more tractable 
than the general two-body problem in general relativity. 

This limit is also of great astrophysical interest, as it 
perfectly describes capture binaries: binary systems cre- 
ated by the capture of stellar mass compact objects onto 
relativistic orbits of massive black holes in galaxy cores. 
Post formation, the evolution of such binaries is driven 
by gravitational-wave (GW) emission — the GW back- 
reaction circularizes and shrinks the binaries, eventually 
driving the smaller body to plunge and merge with its 
larger companion. Such events are now believed to be 
relatively abundant (see Ref. jl| for up-to-date discus- 
sion and review of the relevant literature) . Since the last 
year or so of the inspiral is likely to generate GWs that 
lie in the low-frequency band of space-based GW anten- 
nae such as LISA 2j, extreme mass ratio inspirals (or 
EMRIs) are key targets for future GW observations. 

This paper is the first in a series whose aim is to de- 
velop adiabatic EMRI waveforms. "Adiabatic" refers to 
the fact that they are computed using an approximation 
to the true equations of motion that takes advantage of 
the nearly periodic nature of the smaller body's motion 
on "short" timescales. This approximation fails to cap- 
ture certain important aspects of the binary's evolution. 



In particular, adiabatic waveforms only incorporate dis- 
sipative effects of the small body's perturbation — effects 
which cause radiation of energy and angular momentum 
to distant observers and down the hole, driving the orbit 
to decay. Conservative effects — effects which conserve 
energy and angular momentum, but push the orbit away 
from the geodesic trajectory of the background spacetime 
— are missed in this approach. It has been convincingly 
demonstrated Q that conservative effects change orbital 
phasing in a way that could be observationally signifi- 
cant. The dissipative-only adiabatic approach to EMRI 
waveform generation is thus, by construction, somewhat 
deficient. 

In our view, this deficiency is outweighed by the fact 
that it will produce waveforms that capture the spectral 
features of true waveforms — a complicated shape "col- 
ored" by the three fundamental orbital frequencies and 
their harmonics. Also, the adiabatic approach is likely to 
produce these waveforms on a relatively short timescale. 
Though not perfectly accurate, adiabatic waveforms will 
be an invaluable tool in the short term for workers devel- 
oping a data analysis architecture for measuring EMRI 
events. In the long term, these waveforms may even be 
accurate enough to serve as "detection templates" for 
EMRI events. Measuring the characteristics of EMRI 
sources will require matching data with as accurate a 
model as can be made, and over as long a timespan as 
possible — perhaps a year or more. By contrast, de- 
tecting EMRI events does not require matching a signal 
with a template for such a long time For the short 
integration times needed for detection, work in progress 
indicates that conservative effects do not shift the phase 
so badly that the signal fails to match a template. What 
shift does accumulate due to conservative effects can be 
accomodated by systematic errors in source parameters, 
allowing detection to occur. (This is discussed in Ap- 
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lendix A of Ref. [1], Chapter 4 of Ref. and Refs. 
^, currently in preparation.) 



B. Our approach to adiabatic inspiral 

The approach which we advocate for building adiabatic 
waveforms uses a hybrid of frequency-domain and time- 
domain perturbation theory techniques. These two tech- 
niques have complementary strengths and weaknesses; by 
combining the best features of both toolsets, we hope to 
make waveforms that are as accurate as possible. Though 
a diversion from the main topic of this paper, this ap- 
proach is a key motivation for our work. We thus ask the 
reader to indulge us as we briefly describe our rationale. 

In the adiabatic limit and neglecting conservative ef- 
fects, the separation of timescales means that orbits 
are, to high accuracy, simply geodesic trajectories of the 
spacetime on short timescales. The orbital decay that 
is driven by backreaction amounts to the system evolv- 
ing from one geodesic orbit to another. Computing the 
effect of radiation reaction thus amounts to computing 
the sequence of orbits through which the system passes 
en route to the final plunge of the smaller body into the 
large black hole [9]. 

A geodesic orbit is characterized (up to initial condi- 
tions) by three constants: energy E; axial angular mo- 
mentum L^; and "Carter constant" Q (see, e.g., jTo| . 
Chap. 33). Computing this sequence of orbits is equiv- 
alent to computing the rate at which these constants 
change due to radiative backreaction. In this picture, it is 
useful to regard each orbit {E, Lz, Q) as a point in an or- 
bital phase space, and to regard the rates at which they 
evolve, {E, Lz,Q), as defining a tangent vector to the 
trajectory an evolving system traces through this phase 
space. Adiabatic radiation reaction thus amounts to cal- 
culating this tangent at all orbits. 

In the extreme mass ratio limit, the smaller body 
moves very slowly through orbit space — it spends many 
orbits in the vicinity of each {E, Lz,Q)- This slow evo- 
lution means that the tangent vector is most accurately 
represented by the average rate at which these constants 
evolve: {{E), {Lz), (Q)), where the angle brackets denote 
an appropriate averaging with respect to the orbits. Such 
an averaging is defined in Ref. [3]. 

Once adiabatic radiation reaction data has been found 
for all orbits, it is straightforward to choose initial con- 
ditions and compute the worldline z{t) which an inspi- 
ralling body follows. In this framework, it is just a 
geodesic worldline with the constants slowly evolving: 

z{t) = Zg,od[E{t),Lz{t),Q{t)] . (1.1) 

This worldline can then be used to build the source term 
for the wave equation, allowing us to compute the gravi- 
tational waves generated as the small body spirals in. We 
note here that this approach is conceptually identical to 
the "kludge" presented in Ref. [11]. Indeed, the almost 



unreasonable success of kludge waveforms served as an 
inspiration for this formulation of inspiral^ . 

In the hybrid approach, a frequency-domain code 
would be used for the adiabatic radiation reaction, and 
a time-domain code used to generate the waves from a 
small body following the worldline that radiation reac- 
tion defines. Since any function built from bound Kerr 
black hole orbits has a spectrum that is fully described 
by three easily computed frequencies and their harmon- 
ics [H, ! the averaging needed in this prescription is 
extremely fast and easy to compute in the frequency do- 
main. Many harmonics may be needed, but each har- 
monic is independent of all others. Frequency-domain 
codes are thus easily parallelized and the calculation can 
be done very rapidly. In the time domain, averaging is 
much more cumbersome — a geodesic orbit and the ra- 
diation it generates must be followed over many orbits 
to insure that all beatings between different harmonics 
have been sampled. Convergence to the true average for 
quantities like {E) will be slow for generic (inclined and 
eccentric) Kerr black hole orbits (the most interesting 
case, astrophysically) . 

By constrast, building the associated gravitational 
waveform with a frequency-domain code is rather cum- 
bersome. One must build the Fourier expansion of the 
waves from many coefficients, and accurately sum them 
to produce the wave at any moment of time. The benefit 
of each harmonic being independent of all others is lost. 
In the time domain, building the waveform is automatic 
— modulo two time derivatives, the waveform is the ob- 
servable that the code produces. Given a worldline, it 
is straightforward to build a source for the time-domain 
wave equation; one then cannot help but compute the 
waveform that source generates. 

We are thus confident that by using both frequency 
and time-domain perturbation techniques, we can get the 
best of both worlds — letting each technique's comple- 
mentary strengths shine to build EMRI waveforms that 
are as accurate as possible, in the context of the adiabatic 
approximation. 



C. This paper 

Key to the success of the hybrid approach is the de- 
velopment of fast, accurate codes for both frequency and 
time-domain approaches to black hole perturbation the- 
ory. First results from a frequency-domain code which 
can handle generic orbits have recently been presented 
[13], and the last major formal step (understanding the 



The major difference between the hybrid inspiral described here 
and the kludge is that the hybrid inspiral aims to correctly solve 
a wave equation at all points along the orbit. The kludge instead 
uses a physically motivated approximate wave formula based on 
variation of the source's multipole moments, defined in a partic- 
ular coordinate system. 
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adiabatic evolution of Carter's constant Q due to GW 
emission) is essentially in hand The frequency- 

domain side of this program is thus in a good state. Our 
goal now is to develop time-domain tools sufficiently ro- 
bust and generic to handle the case of interest. 

The major difficulty in building a time-domain pertur- 
bation code is the source term, representing the smaller 
member of the binary which perturbs the large black 
hole's spacetime. In the frequency domain, the small 
body is usually approximated as having zero spatial ex- 
tent, and can be represented using delta functions (and 
their derivatives). One then constructs a Green's func- 
tion from solutions of the source-free perturbation equa- 
tion and integrates over the source. Thanks to the delta 
nature of the source in this representation, this integral 
can be done analytically. This trick cannot be done in 
the time domain — one must choose a functional form 
of the source which can be represented on a finite differ- 
ence grid. The challenge is to pick a representation that 
accurately captures the very narrow spatial extent of the 
source, but is sufficiently smooth that the source does 
not seed excessive amounts of numerical error. This is 
particularly difficult for sources representing highly dy- 
namic, generic Kerr black hole orbits in which the source 
rapidly moves across the grid. 

Much recent success in this approach has come from 
representing the source as a truncated, narrow Gaussian 
Khanna [11] and Burko and Khanna [l^ have so 
far examined some orbit classes (equatorial orbits, both 
circular and eccentric) and found that they can quickly 
and robustly generate waveforms from orbits around Kerr 
black holes. As a diagnostic of this technique, they com- 
pute the flux of energy carried by this radiation and find 
agreement with pre-existing frequency-domain calcula- 
tions at the few percent level. 

An interesting recent development is the use of finite el- 
ement techniques to represent time-domain sources. Such 
methods are tailor made for resolving problems with 
multiple lengthscales, and as such may be ideal for the 
EMRI problem. Sopuerta and Laguna [l^] have found 
that a finite element code makes it possible to repre- 
sent the source with amazing accuracy — agreement with 
frequency-domain calculations at the few hundredths of 
a percent level seems common. To date, they have only 
examined binaries in which the larger black hole is non- 
rotating, but they argue convincingly that the diffi- 
culties required to model Kerr perturbations should not 
be terribly difficult to surmount. These techniques are 
an extremely promising direction that is sure to develop 
extensively in the next several years. 

Our goal in this paper is to develop another represen- 
tation of the source term that is simpler (and concomi- 
tantly less accurate) than finite element methods, but 
that is developed somewhat more systematically than the 
truncated Gaussian. The key ingredient of this approach 
is an extension of a finite impulse representation of the 
Dirac delta function [H, [gj- In essence, one writes the 
discrete delta as a series of spikes on the finite differ- 



ence grid, with the largest spike centered at the argu- 
ment of the delta, and with the spikes rapidly falling off 
away from this center. One chooses the magnitude of the 
spikes such that the delta function's integral properties 
are preserved, particularly the rule that 

J dxf{x)Six-xo) = f{xo) . (1.2) 

The discrete delta described in Ref. 23] allows one to 
make a tradeoff between localization and smoothness — 
one can smear the delta over k points, choosing k to be 
small if source sharpness is the key property needed, or 
allowing k to expand if too much sharpness causes numer- 
ical problems. This representation introduces a kind of 
optimization parameter which one can engineer as needed 
to find the best compromise between smoothness and lo- 
calization. 

We extend the finite impulse representation of the 
delta described in [2^ [2^ in two important ways. First, 
the source term of the Teukolsky equation requires not 
just the delta, but also the delta's first and second deriva- 
tives. We therefore generalize this procedure to develop 
discrete delta derivatives. If the delta is represented by 
k points, then both derivatives will require k + 2 points. 
The guiding principle of this extension is again the no- 
tion that the integral properties of these functions must 
be preserved: 

J dx f{x)5'{x - xo) = -/'(xo) , 

J dxf{x)S"{x~xo) = f"{xo) . (1.3) 

(Here, prime means d/ dx.) 

If the discrete delta function does not lie precisely on a 
grid point, then one must use interpolation to appropri- 
ately weight impulse functions from the neighboring grid 
points. Our second extension of Ref. [23| is to introduce 
higher order interpolation (cubic) which offers another 
way to trade smoothness for localization. This is partic- 
ularly valuable when (as in our application) the source is 
coupled to a wave equation. 

We test this representation by developing a new time- 
domain Teukolsky equation solver which uses this form of 
the delta for its source (the "5-code" ) and comparing to 
a well-established code (see, e.g., [I3,[ll,[lll) which uses 
a truncated Gaussian (the "G-code"). The G-code has 
been described in detail in a previous publication p^j : 
for the purpose of this paper, the most salient feature 
of this code is how it represents the source term. The 
G-code begins with the following approximation to the 
Dirac delta function: 

,[,_,„,._J^„p(_tfW)^ (,4, 

[Cf. Ref. Tt^, Eq. (19).] The width a is chosen to be 
small enough that this delta only spreads across a few 
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grid zones. The Teukolsky equation source is then built 
from this Gaussian representation and its derivatives. 

The (5-code by contrast uses the representation de- 
scribed in detail in the following sections of this paper — 
a representation that is discrete by design, rather than a 
discretization of a continuous delta approximation. The 
principle advantage of this form seems to be that it makes 
it possible to rigorously enforce integral identities involv- 
ing the delta plus its derivatives. 

There are a few other minor differences between these 
two codes, which are artifacts of the codes' independent 
developments. Chief among these differences are the use 
of slightly different axial coordinates (the G-code uses the 
usual Boyer-Lindquist coordinate (f); the ^-code follows 
Ref. and uses a coordinate (f> defined in Sec. Ill A[) . 
and the use of slightly different fundamental "fields" (i.e., 
slightly different representations of the Weyl curvature 
scalar ip4 which the Teukolsky equation governs). There 
are also some differences in the way the two codes imple- 
ment boundary conditions. We present a detailed com- 
parison of the results from the two codes in Sec. IIVI It's 
worth pointing out that that we also have taken the G- 
code and replaced its source term with that used by the 
i5-code. This exercise confirmed all of the results we ob- 
tained with the (5-code, demonstrating that these minor 
differences had no impact on our results. 

As a proof-of-principle check of this idea's validity, we 
restrict our present analysis to circular, equatorial or- 
bits. The results from both codes are then compared 
against frequency-domain results. Flux of energy car- 
ried by gravitational waves is a very useful benchmark 
with which to diagnose a perturbation theory code's ac- 
curacy (especially for very simple orbits when averaging 
is easy both in time and frequency domains). In all cases, 
we find (after some experimentation to optimize our dis- 
crete delta) that this new source form is more accurate 
(typically by factors of 2 — 5) and faster (often by factors 
of about 10) than the truncated Gaussian. For our pur- 
pose, it appears that this form of source function will be 
very well-suited to serve as the core of the time-domain 
portion of our hybrid approach to EMRI waveforms. 

Future papers in this series will then apply this tech- 
nique to flesh out the hybrid approach. Our first fol- 
lowup will examine how well this source works for highly 
dynamical trajectories — generic (inclined and eccentric) 
geodesic orbits, plunging orbits, and non-geodesic trajec- 
tories (standing in for orbits that evolve due to radiation 
reaction). Early results from this analysis indicate that 
the discrete delta source term handles such orbits very ro- 
bustly, validating earlier results for eccentric equatorial 
orbits [.19j; work is in progress to extend this to generic 
orbits. We will then begin developing hybrid EMRI wave- 
forms in earnest, using frequency-domain tools to com- 
pute the effects of radiation reaction, building an inspiral 
worldline from those effects, and finally computing the 
waveform with our time-domain code. 



D. Organization of this paper 

The remainder of this paper is organized as follows. 
Section |ll] reviews how one solves the Teukolsky equa- 
tion in the time domain, introducing the equation itself, 
specializing to the form that we use for our calculations, 
and showing how to extract waveforms and fluxes from 
its solutions. We first review in Sec. Ill A] how one solves 
for the homogeneous (source-free) form of the Teukolsky 
equation, an important first step to developing a robust 
solver for the sourced case. We follow very closely the 
procedure laid out in Ref. [1^ ; this section is thus largely 
a review and summary of that paper (with a few minor 
corrections noted). Section fll Bl then describes in detail 
the form of the source term that applies when perturba- 
tions arise from an orbiting body. 

The need to model this source using a delta function 
motivates Sec. IIIII our model for a discrete delta and its 
derivatives. This section presents the key new idea of 
this paper. After describing the basic idea behind our 
discrete delta, we first present in some detail (Sec. IIII A[) 
an extremely simple two-point discrete delta function. 
This illustrates the concepts and principles of this ap- 
proach. We then generalize this idea to a multiple point 
delta in Sec. IIIIBI and then show how to smooth things 
with higher order interpolation in Sec. IIII CI Some pre- 
liminary issues related to the convergence of quantities 
computed using the discrete delta are introduced in Sec. 

We test this delta representation in Sec. IIVI exam- 
ining how the various methods we develop work at de- 
scribing the Teukolsky source function. Section IIV Al 
first compares the different discrete delta functions with 
each other, demonstrating how the different approaches 
change the quality and accuracy of our results. Based 
on this analysis, we choose to use the high order (cubic) 
delta described in Sec. IIII Cl in the remainder of our work. 
We then examine the convergence of our code, demon- 
strating second-order convergence in Sec. IIVBI Finally, 
in Sec. IIV Cl we compare the discrete delta with the Gaus- 
sian source function, demonstrating explicitly how this 
new representation improves both the code's speed and 
accuracy. 

Our benchmark for evaluating our results is to com- 
pare the energy flux carried by the system's emitted 
gravitational waves to results obtained using a frequency- 
domain code [U, [l^l . This operation requires us to ex- 
tract these waves at a particular finite radius. Section 
IVl examines the dependence of these fluxes as a function 
of extraction radius, and finds that they are very well fit 
by a simple power law. Using this law, we can easily ex- 
trapolate our results to very large radius; doing so greatly 
improves agreement with frequency-domain results, typ- 
ically indicating that our errors are significantly smaller 
than 1% for a large fraction of parameter space. 

Concluding discussion in given in Sec. IVII Besides 
summarizing the major findings of this analysis, we dis- 
cuss in some detail future projects to which we intend to 



5 



apply this new computational technology. 



II. NUMERICAL IMPLEMENTATION OF THE 
TEUKOLSKY EQUATION IN THE TIME 
DOMAIN 



The 9 and (j) directions are taken with respect to the black 
hole's spin axis; the function $(t, r, 9) is a reweighting of 
the field \1/ which we define precisely in Eq. (|2.8p below. 



A. Homogeneous Teukolsky equation 



Here we describe the evolution algorithm used in the 
(5-Code, built using a two step Lax-Wendroff algorithm. 
Our notation and approach closely follow that used in 
[1^; some of this section therefore can be considered a 
summary of that paper. All details related to the G-Code 
were described in p7l |. 

Teukolsky derived a master equation that describes 
perturbations due to scalar, vector and tensor fields in 
the vicinity of Kerr black holes in [13, H^. In Boyer- 
Lindquist coordinates, this equation reads 



Reference [25| demonstrated stable numerical evolu- 
tion of (|2.1[) for s — —2. The (5-code has been built 
using the algorithm presented in [25j . after accounting 
for some typographical errors, which are also discussed 
in [s^l- The contents of this section are largely review 
of the results presented in [1^ ; as such, our discussion is 
particularly brief here. 

Our code uses the tortoise coordinate r* in the radial 
direction, and azimuthal coordinate 4>\ these coordinates 
are related to the usual Boyer Lindquist quantities by 
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where M is the mass of the black 
lar momentum per unit mass, A 

. = M ± V — and s is the 
"spin weight" of the field. The s = ±2 versions of these 
equations describe perturbations to the Weyl curvature 
tensor, in particular the radiative degrees of freedom ■00 
and -04. That is, 'J = for s = +2, and 'i' = p^'^i'i for 
s = -2, with p = -l/(r - mcos6'). The T in the RHS 
of this equation depends on the details of the perturbing 
source. It is here that the Dirac delta function and its 
derivatives enter. A discussion of T is postponed to the 
latter half of this section, after we discuss the numerical 
evolution of the homogeneous Teukolsky equation. 

Gravitational waves, h+ and hx as well as the energy 
flux dE/dt [H, [11], can be obtained far away from the 
system by using s = — 2 in Eq. (|2.ip and then identifying 
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Following [25|, we factor out the azimuthal dependence 
and use the ansatz. 
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allows the Teukolsky equation to be rewritten as 
dtu -f Mdr*u + Lu + Au = T, 

where 

u = {$fl,$/,nfl,n/} 
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is the solution vector. The subscripts R and / refer to 
the real and imaginary parts respectively. (Note that 
the Teukolsky function ^' is a complex quantity.) The 
matrices iVf, A and L are 



M 



b 


"131 





b 

m-32 



-77132 7)7.31 



0\ 

-b 
-bj 



(2.14) 



6 



/ 


0-31 0,32 

-0-32 flSl 




(2.15) 



-a34 033 



and 




where 



"I31 
TO32 

031 



032 

033 
034 



~bci + bdr'b- 
bcj, + 2am{r^ 



C2 



a2)/E2 



(2.16) 



(2.17) 
(2.18) 



= A- 



— sm 



-6A 



Y? sin^ e 

a2 + r (r (s + 2) - M {s + 3)) 



AM{r ~ \)smaM + (6amA)/r 



ci , 
-C3 , 
_ A_^ 

2 I 71,f„2 I „3 I „„2^ 



cot 6* 



A d 



2s{-iMr^ + Ma'' + r'^ + ra^)/I]^ 



^31 

ci 

C2 

C3 = 2a(2rMm + Ascos^)/^ 



= -2' 



^A(l 



- (a^ - r2)A/s 6A6 



£2 



(2.19) 

(2.20) 

(2.21) 
(2.22) 

(2.23) 

(2.24) 

(2.25) 

(2.26) 



The equations above have been written such that the 
typographical errors in 2§|'s 031,032,034 and C2 are ob- 
vious. It turns out that the coefRcients hsted in [2^ are 

correct when the ansatz ^{t,r* ,9,4>) = e''"'^$(t, r*, 6^) is 
used. T is a quantity contructed from the source term 
and is discussed in the latter half of this section. 
Rewriting Eq. (P?T^ as 
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and subjecting it to Lax-WendrofF iterations produces 
stable time-evolutions. Each Lax-WendrofF iteration con- 
sists of two steps. In the first step, the solution vector 
between grid points is obtained from 
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This is used to compute the solution vector at the next 
time step. 
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The angular subscripts are dropped here for clarity. All 
angular derivatives were computed using second order 
centered finite difference expressions. Notice that the 
matrices £), A and M are time independent. In addition, 
the time stepping must satisfy the Courant-Friedrichs- 
Lewy condition [2^, 5t < ma,x{ Sr*, 5M 59}, where St is 
the time step^. 

Following [2^ , we set $ and 11 to zero on the inner and 
outer radial boundaries. While the aymptotic behavior 
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makes this condition reasonably accurate at the inner 
boundary, it is clearly unphysical at the outer boundary. 
By placing our outer boundary sufficiently far, error due 
to our outer boundary condition can be made unimpor- 
tant; reflections from the outer boundary have no impor- 
tant impact on our results. Symmetry of the spheroidal 
harmonics is used to determine the angular boundary 
conditions. For even \m\ modes, we have dg^ = at 
9 ~ O.TT. On the other hand, $ = at 6* = 0, tt for modes 
of odd |m|. 

As a test of our evolution equation, we have examined 
source-free field evolution (setting T = 0) for a variety of 
initial data, in particular comparing extensively with the 
results of [2^. As an example of one of our tests, we take 
initial data corresponding to an ingoing, narrow Gaus- 
sian pulse. This data perturbs the black hole, causing it 
to ring down according to its characteristic quasi-normal 
frequencies. We find extremely good agreement in mode 
amplitude and evolution (typically ~ 1% error) with re- 
sults from ^25*1. Figure [T] shows the result of such a test, 
illustrating the quasi-normal ringing and power law tail 
for the Z = 2, m = mode of a black hole with spin 
parameter a = 0.9. 



B. The source term 

We now consider the source term, T, of Eq. (2.1). It 
is given by 

T = 2p-'*T4 , (2.33) 
Ti = (A + 37-7 + 4At + /i)(A + 27-27 + ^)r™^ 
~iA + 3-f~j + A^ + fl){S-2f + 2a)T„^ 
-t-(^ - f + /3 + 3a + 47r)(^ - f + 2/3 -I- 2a)r„„ 
-(^ - f + ^ + 47r)(A + 27 + 2/2)r„™ . (2.34) 



^ Conducting a von Neumann local stability analysis on all the 
points of our numerical grid yields that this condition is sufficient 
for stable evolutions. See reference |25|| for more details. 
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Quasinormal decay ~ exp(-0.078193 t/M). 




200 250 

time/M 



450 



Fourier transform of tfie decay, normalized sucfi tfiat Maxfl>(f)]=1 . 
Center frequency is (0= 0.41 41 7 IVI"^ . 




2 2.5 3 

angular frequency, coM 



FIG. 1: Illustrations of quasi-normal ringing for a black hole with a/M = 0.9; the / = 2, m = mode is shown here. Top 
panel: Evolution of the magnitude of the Teukolsky function, extracted at r = 20M, 9 — Tr/2. We plot the time evolution 
of In |$| at this position. Overplotted on this curve (dashed line) is a function oc exp(— O.Q78193t/A/), demonstrating that we 
recover the expected decay law with a damping time r = 12.789M. Bottom panel: Magnitude of the Fourier transform of 
Notice that it peaks at a; = 0.41417/M. These results for uj and t are in excellent agreement with the expected values 
of [ujM, M/t) = (0.41, —0.078) from Ref. [2^ for quasi-normal ringing of the I = 2, m — Q mode for a = 0.9. 



Reference [22| provides definitions for the various quan- 
tities which appear in Eqs. (|2.34p and (|2.34p . Of partic- 
ular importance are the quantities T„„, Tnm, and Tmm, 
given by contracting the stress energy tensor for the or- 
biting body with the Newman-Pcnrosc null-tetrad legs 
and m^: 



2 ' (2-^^) 



\/2{r — ia cos 9) 
{—iasm9, 0, 1, —i csc9) 



(2.36) 



where A = — 2M r + and p ^ = + cos^ ' 



operators A and 6: 
d 



A = 



aimp 



dt 



2 dr 



2 



dx^^ 

ia sm9 {r + ia cos 9) d 
di 

(r + ia cos 9) p'^ d mp'^{r + ia cost 



(2.37) 



_ _ ^ 

St+(>8 + 



d9 



V2sm9 



(2.38) 



(The operator A is normally written without the tilde; 
we have added it here to avoid confusion with A = — 
2Mr + a^.) 

To proceed, we next must analyze the stress energy 
tensor describing the small body. A point body of mass 
p disturbing the Kerr spacetime is given by 



Also of great importance here are the Newman-Penrose 



T^. = ^"^"'^^ S[r ~ R{t)]S[9 ~ e{t)]S[ct> - $(t)] , (2.39) 



8 



where it'^ = dx^ /dr, and where B(i), $(<)] describe 

the Boyer-Lindquist coordinate worldhnc of the small 
body. Due to axial symmetry of the Kerr spacetime, 
the only (j) dependence in the stress-energy tensor comes 
from 5[<p ~ ^{t)]- Writing 



im(j} 



(2.40) 



and using the fact that 



5[4> - $(<)] = ^Y^ e™l'^-*(*)l , (2.41) 



we find 



^ _ R{t)]5[e - e(i)]e-™*(*) . (2.42) 

sm 

Using this expansion for T^^, it is a simple matter 



to construct r„„ = n^^n^Tfj^^, T^m 



n^m'^T^^, and 



Tfiifn — 'fn^^fhyT^^. We then insert these terms into Eqs. 
(|2.34p . Using the chain rule repeatedly leaves us with 
a (rather complicated) expression involving radial and 9 
derivatives of Dirac delta functions. We thus face the 
task of representing the delta function and its deriva- 
tives accurately on a numerical grid. This is the major 
innovation of this paper, and is discussed in detail in the 
following section. 

It should be noted at this point that, since the Teukol- 
sky equation is most naturally written in terms of the 
tortoise coordinate r* , we must describe the radial be- 
havior of the source term in r* as well. To this end, we 
replace the radial delta function and all radial derivatives 
as follows: 



5[t - R{t)] 

d_ 

dr 



5[r* -R*{t)] 
\dr/dr*\ 

r"^ +0^ d 
A dr* 



(2.43) 
(2.44) 



III. THE DISCRETE DELTA FUNCTION AND 
ITS DERIVITATIVES 



As pointed out in the previous section, the Dirac delta 
function enters the Teukolsky equation because we ap- 
proximate the perturbing mass by a point particle. By 
its definition as an integrable singularity, the delta func- 
tion is very difficult to represent on a finite difference 
grid. The best we can hope to do is to develop a model 
function that captures its most important features, par- 
ticularly localization to a very small spatial region, as 
well as integr ability and derivative properties. The fol- 
lowing three subsections describe the model for the delta 
function we have developed. We first describe a very ba- 
sic model that demonstrates how to satisfy our criteria in 
SecHEl In Sec. iniB] and Unci we then refine this ba- 
sic model. These refinements have been found to improve 
the overall accuracy of the code. 

The discrete delta function approach we use is inspired 
by the work presented in Rcf. [23] . Our technique can be 
considered an extension of that used in (23j ; in particu- 
lar, they do not develop delta function derivatives, nor do 
they implement all the refinements discussed in llll Bl and 
nil CI Nonetheless, Ref. [23| played an extremely impor- 
tant role in developing the foundations of our work. 

Before turning to a detailed discussion of our tech- 
niques for modeling the delta function on a numerical 
grid, we first mention some general considerations per- 
taining to delta functions on a finite difference grid. For 
concreteness, consider a function and delta combination, 
f{x)S(x — a). The function f{x) is taken to be known, 
and can be calculated for any x. For the sake of argu- 
ment, let the delta be modelled by two coefficients Sk and 
(5fc+i on grid points Xk and Xk+i respectively; the delta is 
taken to be zero everywhere else. (This in fact pertains 
to the form of the delta discussed in Sec. IIII Al ) 

Now imagine integrating f{x)5(x — a) over all x. An- 
alytically, we know that this should give us /(a). Nu- 
merically integrating this on our grid gives us 



Finally, in our numerical implementation, we define 
the vector T appearing in Eq. (2.7) as 



0,0,Re(r),Im(T) 



(2.45) 



h'^J{xi)5{xi- a) = h[f{xk)Sk 

i 

+ /(.Tfc+l)4+l]. (3.1) 



where 

f : 



47rA {r 



exp 



■In 



r — r_|_ 



(2.46) 



The exponential factor in this expression corrects for the 
fact that the evolution code uses the azimuthal variable 
(j), but the source term is expanded in (j). 



This equation suggests that the numerical integral ap- 
proximates f{a) by interpolating between grid points Xk 
and Xk+i- If f{x) is rapidly varying, this interpolation 
may not be accurate enough; this is sure to be a source 
of error as we integrate our PDE forward in time. 

Great improvement can be achieved by enforcing the 
well-known identity 



fix)S{x - a) = f{a)5{x - a). 



(3.2) 
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The numerical integral now becomes 

h^f{a)5{xi-a) = h[f{a)Sk + f{a)dk+i] 

i 

= f{a)h{Sk + Sk+i) 

= /(a) , (3.3) 

and the identity is preserved exactly. In the last step, we 
use the discrete analog of the property 



/ 



dx S(x — a) = 1 



(3.4) 



Similar identities can be used on the delta function 
derivatives: 

f{x)d'{x-a) = f{a)S'{x~a) 

-f{a)S{x-a), (3.5) 
f{x)S"{x-a) = fia)6"ix-a)~2f'{a)S'{x-a) 

+f"{a)6{x^a) . (3.6) 

We recommend using these identities as much as 
possible when numerically implementing the algorithms 
sketched in the following three subsections. 



A. A simple numerical delta function 

Consider the function S{x — a), where Xk < a < Xk+i; 
i.e, a lies between two discrete grid points. Let h — 
Xk+i ~ Xk — Xk — Xk-i be the grid resolution. We use 
the following integral to define the delta function: 



dx f{x) 5{x — a) = f(a) , (3-7) 



where e > and / (x) is any well behaved function. This 
means that 5{x — a) is zero everywhere, except at x = 
a, where it is singular. Translating this integral to a 
summation, we have: 



f ^ dxf{x)6{x-a) ~ hY,fi^^)^^ (3.8) 

Ja-e j 



where Si is the discrete delta defined on the grid. Since 
a does not necessarily lie on a gridpoint, we can linearly 
interpolate to find: 



r( ^ f{xk+l) - f{xk), \ , fl \ 

/(a) = {a-Xk) + f{xk) 

+0{h^) . (3.10) 



Substituting this back into our earlier expression and 
comparing coefficients, we have 



5^ = 



a - Xk 

Xk+l - 



for i = k 
— tor I — 



1 



everywhere else 



(3.11) 



Notice that if a = Xk, then Si = l/h for i = k, but is zero 
everywhere else; a similar result holds if a = Xk+i- This 
reproduces our intuitive notion that the delta function 
is zero everywhere except at a single point, and that it 
integrates to unity. We take the viewpoint that the inte- 
grability of the delta function is its key defining property, 
using this rule to derive the results presented below. This 
is the approach that was used in Ref. [23j . 

Another approach to defining a numerical delta func- 
tion, suggested in [l^, is to first define a step function 
on the grid, and then use finite differencing to obtain the 
delta and its derivatives. The approach described above 
matches this proposal when a lies exactly on a grid point. 

We can proceed in a similar fashion to find formulae 
for the derivatives. Let us define 
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Xk+i - a 



1 



a - Xk 



h h 
Again, we start from the defining integrals, 

dx f (x) S' {x ~ a) ~ ~f'{a) 



(3.12) 



(3.13) 



j dx f{x)S"{x-a) = /"(a). (3.14) 

Note that a prime denotes d/dx. Our goal is to derive a 
form which enforces these integrals in summation form: 

fix,) SI c ~f'{a) 

i 

= -hY,i'{x.,)s, 

i 

= -if'ixk) 

-(l-7)/'(^fe+i) + 0(/»');(3.15) 

hj2 /(2^.)^r ^ /"(«) 

i 

i 

= lf"{xk) + 

{\~^)f"{xk+i) + 0{h^).{?..lQ) 

The derivatives of f{xk) are given by the finite difference 
formulae, 

f (xk) = + 0{h), (3.17) 



2h 

ri, I X f {xk+i) - 2f {Xk) + f {Xk-l) 

+0{h^). 



(3.18) 
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Substitution of these approximations in (|3.15p and (|3.16D 
and a comparison of coefficients yields for the derivative: 



J_ 

2h' 

1-7 

!_ 

7-1 



for i = k , 
for i — k + 1 
for i ~ k + 2 



everywhere else 



(3.19) 



For the second derivative: 



S'/ = -I for i = fc - 1 , 



1 - 37 
37-2 

1-7 



for i — k , 
for i ~ k + I 
for i — fc + 2 , 



/i3 

everywhere else 



(3.20) 



Notice that we need four points to represent the deriva- 
tives of the delta function in this scheme. 



B. A multiple point delta function 



The procedure described in Sec. IIII Al can be extended 
to represent the delta function over a larger number of 
points. On the one hand, this spreads out the delta, 
moving us away from our ideal of a function that is non- 
zero in as small a region as possible; on the other hand, it 
allows us to represent it more smoothly on our grid. The 
number of points (2n -I- 2) that we use can thus be con- 
sidered an optimization parameter, allowing us to trade 
localization for smoothness. As we shall see, there is 
typically a value of n that represents a very good com- 
promise. 

We start off with the 'linear hat' delta function defined 
in [H and [H 



where 



Ji/h for \xi — a\ < e — nh 
otherwise 



7. = -(l-^ 
n \ nh 



(3.21) 



(3.22) 



and n is an integer. Note that when n = 1, 7^ reduces 
to the 7 that was defined in Sec. IIII Al Note also that 
7j is non-zero only for z G [k, ...,k + 2n — 1] (so that there 
are a total of 2n points) , and that a lies between the grid 
points Xk+n-i and Xk+n- In this labeling scheme, Xk is 
the smallest gridpoint where 5i is nonzero. 



Substituting this form of the delta function into our 
defining integral relation. 



we find 



^/(a;,)7.^/(«), 



(3.23) 



(3.24) 



The quantity 7^ is thus a weighting factor whose weight 
depends on the distance of Xi from a. Setting f{x) = 1, 
we find the property 



E7» = l 



(3.25) 



Now consider the derivative of the delta function. Our 
goal is again to enforce the rule 



(3.26) 



Inserting the finite difference formulae for the derivatives 
of /, Eqs. (I3.17P and (|3.18p . into this relation, we find 



/'(«) 



fe+2ri-l 



■i—k 



fjx.+i) - f{x^-.i) 
2h 




1 

2h 



"2h ['>''=+2"-2/(^fe+2n-l)] 

~2h ['>''=+2"-i/(-^fe+2")] 



1 

'2h 

1 



E 7»+i/(a;«) 
^ ■i=fc+i y 



+ ^ [-lkf{xk-i) - 7fc+i/(a;fe) 



Comparing coefficients, we read off 



51 



SI 



k+2n-l 



Si 



k+2n 



Ik 

2/l2 ' 

Ik+l 

2h^ ' 

7fc+3-l - 7fc+3 + l 

2h^ 

for j e [l,2n- 2] 

7fe-|-2n-2 

2h^ ' 

7fe+2n-l 

2h^ ■ 



(3.27) 

(3.28) 
(3.29) 

(3.30) 
(3.31) 

(3.32) 



The formulas for the delta derivative coefficients can 
be understood intuitively. The n-point generalization 
approximates the delta function as an isoceles triangle 
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centered at a and sampled at 2n points. The derivative 
is simply the slope of this isoceles triangle at all points, 
except at the center and the edges, where the derivative 
is discontinuous. The discontinuity is replaced by coef- 
ficients that ensure the integral properties of the deriva- 
tive. The delta derivative takes a particularly simple for- 
mula in the "bulk" : 



Notice that the second derivative is zero in the "bulk" 
— it corresponds to the second derivative of a line, with 
constant slope. Like the first derivative, these coefficients 
are non-zero for i € [fc — 1, fc -|- 2n\ — two points broader 
than the delta itself. 



Si 



k+j 



7fc+j-i 



2nh2 



P:^±^ for je [l,2r 

Xk+j+l " 



2h2 

\xk+j-i - a\ 



2] 



nh 



nh 



for Xk+j+i < a — h 
for Xk+j+i > a + h 



(3.33) 



Notice that the delta derivative coefficients are non-zero 
for i (z [k — l,k + 2n\ — one point wider in each direction 
than the span of the delta on the grid. 

A similar analysis can be done for the second deriva- 
tives. We start off with 



/"(«) 



k+2n-l 

E 

i—k 
A:+2n-l 

E ' 

i—k 



l^f"{x) , (3.34) 
/(.x,+i)-2/(a;,) + /(x,_i) 



We have found that a very sharp delta function, like 
the two-point model described in the previous section, 
leads to instabilities for orbits with varying r or 9 (e.g., 
for eccentric orbits) . Using a smoother n-point represen- 
tation suppresses these instabilities; this will be discussed 
in greater detail in [3ll |. Since r and 6 are constant for 
circular, equatorial orbits, these instabilities do not arise 
for the cases examined in detail here. Thus for the case 
at hand, our numerical errors originating from the ffnite 
representation of the delta are smallest when n = 1. This 
is demonstrated in detail in Sec. IIVI 

C. Higher order interpolation for smoothness 



(3.35) 



Reading off the coefficients leaves us with 



X" 



X" 



X" 



X" 

"k+2n- 



A" 

"k+2n 



7fc4 



Ik+j+i - ^Ik+j + 7fc+j-i 
for j e [l,2n- 2] , 

7fc+2n-2 — 27fe+2ri-l 



/l3 



7fc+2r7 



/l3 



(3.36) 
(3.37) 

(3.38) 
(3.39) 
(3.40) 



Finally, we present a representation of the delta which 
uses a higher order interpolation scheme. This again 
spreads the "stencil" of the delta function over a wider 
patch of the grid, but improves our ability to reproduce 
the integral formulation of the delta identities. 



Using cubic interpolation (which requires a total of 
four points), we find the rule 



i 



/(«) 

(g - Xfc+i)(a - Xfc+2)(a - Xk+3) 

(g - Xk)ia - Xk+2)ia - Xk+3) 
2h^ 

{a - Xk){a - Xk+i){a - Xk+s) 
2h3 

(g ~ Xk)ia - Xk+i)ia - Xk+2) 
6/i3 



(3.41) 



f{xk) 



f{Xk+l) 



+2} 



f{Xk+3.) ■ 



(3.42) 
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The location a lies between grid points Xk+i and Xk+2- From this expression, we read off the coefficients 

, (a-x,+i)(a-^.+2)(a-x.+3) ^^^^^ ^3 43^ 

at Xfc+i , (3.44) 

at Xk+2 , (3.45) 

at Xk+3 ■ (3.46) 







6/l4 




(a 


- Xk)ia 


- Xk+2)ia 


- 2^fc+3) 










(a 


- Xk)ia 


- Xk+i){a 


- Xk+l.) 










{a 


- Xk)ia 


- Xk+i){a 


- Xk+2) 





A similar analysis for the first derivatives yields 



S[ = i->^+i--)i-~^^^+^)i^~^>^+^) , (3.47) 

_ (a- xk){a ~ xk+2){a^xk+,) ^3 



4/i5 

(g - Xk+i){2a - 3xk + Xk+2){a - Xk+3) 

(g - Xk){a - Xk+2){2a + 2:^+1 - ixk+a) 

{a - Xk){a - Xk+i){a - Xk+3) 
4/i5 

(g - Xk){a - Xk+i){a - Xk+2) 



at Xk+i , (3.49) 

at Xk+2 , (3.50) 

at Xk+3 , (3.51) 

at Xk+i ■ (3.52) 



Note, we have used the second order finite difference formula, Eq. (j3.17p . to derive this result. In principle, higher 
order formulas for the derivative could have been used. We kept to the second order formula in order to keep the 
derivative stencil narrow, and also for consistency with our time-stepping algorithm. 

Finally, for the second derivatives we find 

c// jXk+l - «)(« - Xfc+2)(g - Xk+3) . coA 

= ^ at Xk-i , (3.53) 

_ (5g - 3xk - 2xk+i){a - Xk +2){a - Xk+3) . 

at Xk , [o.ui) 



(lOg^ - {9xk + 4.Xfe+i + 7a;fe+2)g + Xk+iXk+2 + 5xk{xk+i + 2xk+2)) {a - Xk+3) 
(g - Xk) (iQg^ - (7xfc+i + 4,Xk+2 + 9a;fc+3)g + 3xk+2Xk+3 + Xk+i{xk+2 + 6xfc+3)) 

(g - Xk)ia - a;fc+i)(5g - 2xk+2 - 3x^+3) 
(g - Xk)ia - Xk+i){a - Xk+2) 



at Xk+i , (3.55) 
at Xk+2 , (3.56) 

at Xk+3 , (3.57) 



^* 2;fc+4 . (3.58) 



As discussed in more detail in the following section, our 
analysis suggests that this cubic interpolation method 
works best. 

We emphasize at this point that, although we are moti- 
vated by Teukolsky equation applications, our discussion 
here was not specialized to the Teukolsky equation in any 
way. The delta models sketched here can be used in any 
finite-difference numerical algorithm. We also note that 
one does not need to stop at cubic-order interpolation; 
the basic idea of that scheme could easily be extended 
to higher order if the application warranted it. As the 



order is increased, the "stencil" of the delta is likewise 
increased, pushing us away from the intuitive notion of a 
structureless impulse. This leads us to believe that there 
may be a certain interpolation order beyond which the 
model ceases to work well. 



D. Convergence with the discrete delta function 

The non-smooth nature of the discrete delta func- 
tion makes understanding the convergence properties of 
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a code based on this function somewhat subtle. Here 
we briefly summarize some key issues related to conver- 
gence with the discrete delta. This summary is based 
on detailed discussion of discretization errors given in 
Ref. 24^. The punchline of this discussion is that the 
discrete delta function is typically at least second-order 
convergent, and thus we expect our code to likewise be 
second-order convergent. 

Let Si be the discretized version of 6{x— a) defined on 
a discrete grid Xi, let h — Xi+i — Xi be the grid spacing, 
and let 5i be non-zero at , x^+i , . . . , Xk+2n-i- The con- 
tinuous variable a lies between Xk+n-i and xu+n- This 
allows us to define a parameter 77 such that 



a = Xk+n-i + vh- 



(3.59) 



This quantity is a measure of how close a is to a grid 
point; clearly, < < 1. 

We now define the moments of the discrete delta by 

fc+2n-l 

Mr{5,a,h) = h 5,{xi-ay , (3.60) 

i—k 

where r is an integer. In the continuum limit, this defi- 
nition becomes 



Mr —> J 5{x — a) (x — a) dx 

= 1 r = 
= r > . 



(3.61) 



A discrete representation will clearly have the correct ze- 
roth moment; however, it will only have Mr>o = up to 
some maximum r. Reference [2^ proves that, if q is the 
lowest integer such that Mq ^ 0, then 



f{a)-hY^5.J{x,) 



< Chi 



(3.62) 



where C is approximately'^ a constant. This delta repre- 
sentation is then gth-order convergent. 

For the multiple point discrete delta discussed in Sec. 
HIlBl we find Mo = 1, Mi = 0, M2 + 0. When we use 
this discrete delta, we therefore expect our code to be 
second-order convergent. We demonstrate this behav- 
ior in Sec. IIVBI For the cubic delta function, we find 
Mq = 1, A/1,2,3 = 0, M4 ^ 0. In this case, errors due to 
the delta representation are expected to be fourth order. 
However, since our stepping algorithm is itself second- 
order, we expect second-order convergence overall. 



IV. NUMERICAL IMPLEMENTATION AND 
EVALUATION OF THE DISCRETE DELTA 
FUNCTION 



We now implement the Teukolsky equation's source 
term using the techniques discussed in Sec. IHI] for the 
simple case of a point particle in a circular, equatorial or- 
bit around a massive black hole. Our goal is to compare 
the different forms of the discrete delta discussed in the 
previous section and to evaluate which is likely to work 
best for practical modeling of radiation from astrophysi- 
cal systems. We also compare our discrete delta model 
to the Gaussian approximation that has been used in 
previous work, illustrating the power of this new model. 

We obviously require some "standard" to compare our 
results against. Frequency-domain codes provide ex- 
tremely accurate results for circular, equatorial orbits, 
largely since their emitted radiation is concentrated in a 
small number of multipoles; as such, they make an excel- 
lent standard against which to compare our results. We 
use the code described in [23] as our standard. 

Tables U and [TTl shows energy fluxes obtained from our 
code for the most dominant azimuthal modes, |r7i| = 2 
and 3 respectively. We compare these figures with those 
obtained from the code used in [27j'^. 

There are two major reasons that the fluxes we com- 
pute depart from those computed by frequency-domain 
codes. First, the time-domain code must extract fluxes 
at some finite radius. The FD approach produces, by 
construction, the waveforms and fiuxes that would be 
measured infinitely far from the generating binary; this 
simply cannot be done on a finite radial grid. A de- 
tailed discussion of the impact of finite extraction radius 
is given in Sec. [V] In brief, we find by varying the ex- 
traction radius that fiuxes can be fit to a very simple 
power law. This power law then allows us to infer the 
fiux that would be measured by distant observers. The 
second source of error is simply numerical — finite differ- 
ence errors plus the approximate nature of our discrete 
delta. Roughly speaking, accounting for finite extraction 
radius reduces our errors by about a factor of 2 - 5; the 
residual error is thus most likely simply numerical error. 
This is described in much greater detail in Sec. [V] 

Figures O and [3] illustrate a typical example of the 
structure for the Teukolsky function that we find. 
We show the to = 2 mode of an orbit with radius 
To — 7.9456M around a Schwarzschild black hole; this 
orbit was selected in order to compare with results pub- 
lished in [lOl . Note that the orbital period at this radius 
is T = 2i^^rllM ~ 140M. The data is read out at 
radius R = Textract = 250M, = 7r/2; our numerical 



^ In our application, C varies slightly depending on how close the 
delta peak is to a grid point. 



* Symmetry in the azimuthal direction results in equal fluxes for 
+m and — m modes. Thus, |m| refers to the sum of the fiuxes 
from the +m and —m modes which is equal to twice the flux 
from either the +m or the —m mode. 
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TABLE L Energy flux extracted at R = rgxtract ~ 250M for circular, equatorial orbits for the m = |2| mode of a particle with 
mass I-l/M = 1. a/M is the BH spin, tq/M is the orbital radius. The labels "5" and "G" refer to the results from 5-code and 
G-code respectively. Values listed under "FD" are the corresponding fluxes from the frequency-domain code used in [27l |. 



a/M 


ro/M 


Es 


Efd 


{Es — Efd) / Efd 


Eg 


{Eg — Efd) / Efd 





6 


7.385 X 10"'' 


7.368 X 10"" 


0.0023 


7.246 X 10"" 


-0.017 





8 


1.650 X 10"" 


1.651 X 10"" 


-0.0055 


1.623 X 10"" 


-0.016 





10 


5.344 X 10"^ 


5.374 X lO""^ 


-0.0004 


5.290 X 10"" 


-0.016 


0.5 


6 


5.551 X 10""* 


5.539 X 10"" 


0.0022 


5.437 X 10"" 


-0.018 


0.5 


8 


1.399 X 10-"* 


1.401 X 10"" 


-0.0015 


1.375 X 10"" 


-0.019 


0.5 


10 


4.781 X lO"'^ 


4.812 X 10"^ 


-0.0065 


4.691 X lO"'^ 


-0.025 


0.9 


4 


2.654 X 10"^ 


2.662 X 10"^ 


-0.0030 


2.611 X 10"^ 


-0.019 


0.9 


6 


4.614 X 10""* 


4.621 X 10"" 


-0.0016 


4.531 X 10"" 


-0.020 


0.9 


8 


1.249 X lO""* 


1.254 X 10"" 


-0.0039 


1.230 X 10"" 


-0.019 


0.9 


10 


4.419 X 10"^ 


4.456 X 10"^ 


-0.0084 


4.339 X 10"^ 


-0.026 


0.99 


4 


2.469 X 10"^ 


2.484 X 10"^ 


-0.0059 


2.434 X 10"^ 


-0.020 


0.99 


6 


4.450 X 10"" 


4.461 X 10"" 


-0.0024 


4.372 X 10~" 


-0.020 


0.99 


8 


1.221 X 10"" 


1.226 X 10"" 


-0.0041 


1.201 X 10"" 


-0.020 


0.99 


10 


4.346 X lO"'^ 


4.386 X 10""^ 


-0.0090 


4.270 X lO"""' 


-0.026 



TABLE IL Energy flux extracted at 250M for circular, equatorial orbits for the m — ]3| mode of a particle with mass /i/M — 1. 
All symbols and notation are as in Table H] 



a/M 


ro/M 


Es 


Efd 


[Es — Efd)/ Efd 


Eg 


{Eg — Efd) / Efd 





6 


1.465 X 10"" 


1.460 X 10"" 


0.0035 


1.431 X 10"" 


-0.020 





8 


2.445 X 10"^ 


2.449 X lO"'^ 


-0.0017 


2.399 X 10"^ 


-0.020 





10 


6.383 X 10"*^ 


6.435 X 10"'' 


-0.0080 


6.291 X 10~" 


-0.022 


0.5 


6 


1.015 X 10"" 


1.014 X 10"" 


0.0011 


0.992 X 10"" 


-0.021 


0.5 


8 


1.993 X 10"^ 


1.980 X lO"'' 


0.0066 


1.935 X 10"^ 


-0.023 


0.5 


10 


5.521 X 10"'' 


5.572 X 10"" 


-0.0090 


5.410 X 10"" 


-0.029 


0.9 


4 


6.485 X 10"" 


6.467 X 10"" 


0.0028 


6.336 X 10"" 


-0.020 


0.9 


6 


8.031 X 10"^ 


8.043 X 10"^ 


-0.0015 


7.865 X 10"^ 


-0.022 


0.9 


8 


1.710 X 10"^ 


1.717 X lO"'' 


-0.0043 


1.677 X 10"^ 


-0.023 


0.9 


10 


4.992 X 10"^ 


5.044 X 10"" 


-0.0103 


4.893 X 10"" 


-0.030 


0.99 


4 


5.932 X 10"" 


5.924 X 10"" 


0.0014 


5.805 X 10"" 


-0.021 


0.99 


6 


7.688 X 10"^ 


7.688 X 10"'' 


-4.8 X 10"^ 


7.528 X 10"^ 


-0.022 


0.99 


8 


1.6542 X lO"'^ 


1.669 X 10"^ 


-0.0086 


1.628 X 10"^ 


-0.025 


0.99 


10 


4.879 X 10"'^ 


4.942 X 10"" 


-0.0127 


4.792 X 10"" 


-0.030 
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grid runs from — lOOAf < r* < 500Af, with a resolution 
Sr* = 0.0625Af, and from 0<9 <n with 69 = n/AO. 

Figure [2] shows the real part of 4" over a broad span of 
time, from the beginning of our simulation to t ~ 800A/. 
At i ~ 250M, a very high amplitude, unphysical burst of 
radiation reaches the extraction radius. This spurious 
burst is due to our initial conditions: We initially set 
^ = and dt^ = 0, which is not consistent with our 
source function. The time at which this burst reaches the 
extraction radius is perfectly consistent with radiation 
propagating at the speed of light (c = 1 in our units) 
across our numerical grid. The burst quickly propagates 
off the grid, and the solution for "if settles down to a 
simple oscillation. This is shown in Fig. [31 which zooms 
in on the behavior of 5' for t > 350M. Notice that -^r 
has an oscillation period of about TOM, precisely what 
we expect for the m = 2 mode of a source whose orbital 
frequency is 140M. The energy flux we find from this 
mode is E/fi^ = 1.708 x 10"*^, in excellent agreement with 
results pubhshed in Ref. [lol (compare Table II of [20j . 
noting that our results require summing over all I for 
fixed m). 




400 500 600 700 800 



Time /M 

FIG. 2: The real part of the m — 2 mode of the Teukolsky 
function as a function of time for a point particle of mass 
fj,/M — 0.01 in a circular orbit of radius tq/M = 7.9456. 
These data were extracted in the equatorial plane [6 — n/2) 
at radius R = 250M. At this location the Teukolsky function 
is zero by construction until t ~ 250M, at which point a 
spurious burst reaches the extraction radius. This burst is 
due to our unphysical initial conditions; it quickly propagates 
off the grid, leaving a reasonable physical solution for all time 
afterwards. 



Figure|3]illustrates the spatial behavior of the real part 
of \[' at a particular moment in time {t = 312M). This 
plot illustrates the behavior of Re ^ as a function r* and 
9 over a wide span of our grid. Clearly visible are the 
m = 2 mode of the radiation propagating to large radius 
as well as the nearly singular delta function source itself. 



Teukolsky function, >? at r =250 and 8 =1 .5708. 




350 400 450 500 550 600 650 700 750 
Time M 



FIG. 3: The same as Fig. [2l but zooming in on the data for 
t > 350 A^. Solid and dashed lines are the real and imaginary 
parts of *!/ respectively. The Teukolsky function oscillates 
with a period of about 70Af; since the source has a period of 
about 140Af, this is exactly what we expect for the m — 2 
mode. We measure the total fiux of energy carried by this 
mode to be E/fi^ — 1.708 x 10~*, in agreement with previous 
results (see, e.g., Ref. j2Q|]). 



A. Comparison of different discrete delta functions 

In this section, we compare results from the various 
models for the delta function presented in Sec. IIIII The 
variable point approximation for the linear delta func- 
tion, presented in Sec. IIIII provides us with a nice handle 
to study the convergence of our results. As we increase 
n, the half-width of the delta function, the singularity 
spreads out and its sharpness decreases. Notice that the 
physical spread of the source term is 2{n + l)Sr*{due to 
the spread of the delta derivatives). Thus, decreasing the 
resolution has the same effect on the physical width as 
decreasing n. 

In Tables Hill and HVl we present results describing two 
equatorial circular orbits, one in the extremely strong 
field {ro/M = 2.32, a/M = 0.9), and another at weaker 
field {ro/M = 12, a/M = 0). We show the variation 
in flux with n, the half-width of the radial delta func- 
tion. The angular delta function is represented using 
two points (i.e., riang = l)j the minimum number of non- 
trivial points. The resolutions (Sr^^S9) are held fixed at 
(0.0625A/, 7r/40) 

The third and fourth columns of Tables IIIII and IIVI 
compare the flux in energy carried by GWs as computed 
using the (5-code to flux computed using our frequency- 
domain standard. The third column gives a "raw" com- 
parison — we extract the time-domain fluxes at radius 
R = 250Af and compare to the frequency-domain result. 
In the fourth column, we extrapolate the time-domain 
data, R ^ oo, using the algorithm described in Sec. |Vl 
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tend to become more important. This can be compen- 
sated by increasing the width of our delta representa- 
tions. At 71 = 7, the spread of our source is just enough 
to accurately compensate for finite difference errors. At 
larger radius (e.g., the vq — 12M orbit shown in Ta- 
ble lIVp , finite difference errors are so small that we do 
best using the minimum number of points possible in our 
model. 



B. Convergence of our code 

As discussed in Sec. IIII D[ we generally expect a code 
built using the discrete delta on grid with spacing h to 
exhibit 0{h?) convergence. In particular, we expect the 
Weyl scalar ■04 to show second-order convergence. We 
check this expectation by examining the flux of energy 
carried by gravitational waves. Since we expect ■04 = 
^true _|_ Q(^fi'^'j^ TffQ likewise expect E to exhibit second- 
order convergence: 



FIG. 4: The same as Fig. 2, but now showing the data for a 
given moment in time {t = 312. 5A/) for a wide range of r* and 
9. Along with the outward propagating radiation packet vis- 
ible at large radius, the nearly singular delta function source 
is clearly visible at the particle's position. 



The fourth column thus contains the most relevant data 
for assessing which delta representation is "best". We 
include the third column to demonstrate that the differ- 
ence before extrapolating is not terribly large, but that 
it is large enough that the gain due to this extrapolation 
is significant. It also illustrates that not performing the 
extrapolation can mislead regarding which form of the 
delta is most accurate. 

Data from Table [TVl indicate that, among the n-point 
representations, n = I gives the best results for circu- 
lar, equatorial orbits. The cubic representation, how- 
ever, is even better — the smoothness of this technique 
apparently reduces error even more. We choose the cubic 
delta for the remainder of our analysis because it is both 
smooth and accurate. 

By contrast, the most accurate flux in Table [1111 occurs 
when n = 7 (total of 8 points to represent the source), 
rather than n = 1. The reason is due to a competition 
between errors from finite differencing and errors from 
our delta representation. In particular, we have noticed 
experimentally that finite difference errors tend, on av- 
erage, to spuriously decrease our measured flux; errors 
from spreading the delta over the grid tend to augment 
the flux. (We emphasize that this is merely a rule-of- 
thumb tendency we have noted; we also emphasize that 
we do not as of yet have a good explanation for these 
effects.) 

As we approach the horizon, finite difference errors 



£;~l04p~K"1' + o(/i') 



(4.1) 



To demonstrate this convergence, we show the energy 
flux measured at i? = 250M for two different strong-field^ 
orbits: ro = 5M, a = 0.8M and ro = 4.64M, a = 0.9M. 
The radial and angular grids are set to 



Sr* = 0.0625 x 2"''/^ , 
S9 = 7r/30 X 2-^/4 



(4.2) 
(4.3) 



Actually, 69 was modified slightly from this to insure that 
6 — tt/2 lies exactly on a grid point. This reduced vari- 
ations about the main trend owing to the slight de- 
pendence of the proportionality "constant" on the delta's 
peak (cf. discussion in Sec. IIII D|) . Figure [5] shows the re- 
sults our runs for 6 G [— 1, . . . , 4]. Convergence is shown 
by examining the fractional error with respect to our 
densest grid, 



\Eb~E4\ 
Ea 



(4.4) 



where _Ef, is the flux inferred at grid parameter h. We 
normalize to 6 = 4 since it is the densest grid available 
to us. Modulo some slight oscillations, the overall trend 
of our data is in very good agreement with second-order 
convergence. 



^ It's worth noting that we found it to be rather difficult to demon- 
strate convergence using weak-field orbits. For such orbits, the 
differences in our computed fluxes were quite small as we var- 
ied our grid density. We need strong-field orbits in order for the 
errors to be large enough that the convergence trend is apparent. 
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TABLE III: Comparison of several implementations of the discrete delta function. We show results for the linear hat delta 
described in Sec. HUB] as well as the cubic delta function described in Sec. lIIICl All fluxes are measured at R — 250M for the 
\m\ — 2 mode. For the linear hat delta, the total number of points in the function is 2(n + 1). The cubic delta uses 6 points 
in all. These results are for orbits of radius rg = 2.32M about a black hole with a — 0.9A/. The total flux in |m| = 2 modes 
according to our frequency-domain standard is Ep-o/jj^ — 2.061 x 10~^. 



Total points, 2(n - 


f 1) 


E250 


(-E250 — Efd) / Efd 


Eao 


{Eaa — Efd) / Efd 


64 




2.889 X 10~ 


2 


4.0 X 10"^ 


2.890 X 10" 


-2 


4.0 X 10"^ 


32 




2.194 X 10" 


2 


6.5 X 10"^ 


2.195 X 10" 


-2 


6.5 X 10-2 


16 




2.055 X 10" 


2 


-2.8 X 10"^ 


2.056 X 10" 


-2 


-2.4 X 10"^ 


8 




2.027 X 10" 


2 


-1.6 X 10"^ 


2.028 X 10" 


-2 


-1.6 X 10-2 


4 




2.023 X 10" 


2 


-1.9 X 10-2 


2.024 X 10" 


-2 


-1.8 X 10-2 


cubic 




2.024 X 10" 


2 


-1.8 X 10^2 


2.024 X 10" 


-2 


-1.8 X 10-2 



TABLE IV: Same as Table Hill but now for an orbit with ro = 12Af about a black hole with a = 0. The frequency-domain flux 
for |m| = 2 modes in this case is Efo/h^ = 2.172 x 10"^. 



Total points, 2(n + 1) 


-E250 


(-E250 — Efd) / Efd 


Eoc 


{Eao — Efd)/ Efd 


64 


2.342 X 10-^ 


7.8 X 10-2 


2.376 X 10-^ 


9.4 X 10-^ 


32 


2.191 X 10-^ 


8.7 X 10"^ 


2.224 X 10-^ 


2.4 X 10-^ 


16 


2.156 X 10-^ 


-7.5 X 10-^ 


2.187 X 10-^ 


7.1 X 10"^ 


8 


2.148 X lO"'^ 


-1.1 X 10-2 


2.179 X lO"'^ 


3.3 X 10-=^ 


4 


2.146 X lO"'^ 


-1.2 X 10-2 


2.177 X lO"'^ 


2.5 X 10-3 


cubic 


2.145 X lO"'^ 


-1.2 X 10-2 


2.177 X lO"'^ 


2.3 X 10-3 



C. Comparison of discrete and Gaussian 
approximations for the numerical delta 

The work in [isL fioj approximates the delta- 
function by a narrow Gaussian such that it integrates 
to unity over the numerical grid. The Gaussian smears 
out the singularity and thus the source term is non zero 
at a large number of points on the numerical grid. In 
contrast, the models presented in Sec. Illll use only a few 
points to depict the delta-function and thus do not share 
this disadvantage. A comparison on the same hardware 
and software platform showed that the techniques used 
here (the 5-code) are about twelve times faster than the 
ones that use a smeared Gaussian (the G-code) . The last 
two columns in Tables U and |TT] show the fluxes from the 
Gaussian-approximated delta function. Note that the er- 
rors in these fluxes are about 2 — 3%, quite a bit higher 
than errors from the (5-code. Both codes were run with 
identical parameters and grid resolutions. The accuracy 
of both codes improves with higher resolution, and im- 
provement in both is seen by moving the extraction ra- 
dius farther out. However, when these parameters are 
fixed, we find that the (5-code is faster and demonstrates 
higher accuracy. 



V. ACCOUNTING FOR FINITE EXTRACTION 
RADIUS 

When one discusses the gravitational-wave fluxes 
which a system generates, one is normally interested in 



their asymptotic value inflnitely far away. It is of course 
not possible for a flnite coordinate grid to reach all the 
way into this distant zone, so it is of great importance to 
understand how our fluxes vary with respect to our finite 
extraction radius R. 

In flat spacetime, the extraction radius is not very im- 
portant; it just needs to be sufficiently far away that the 
field it measures is purely radiative (i.e, not contami- 
nated by near-field effects) . Conservation laws guarantee 
that fluxes follow a 1/r^ law in this region, and so the 
integrated flux is independent of extraction radius. 

Things are not so simple in a curved spacetime — 
radiation is effectively scattered off of spacetime curva- 
ture, modifying its propagation characteristics compared 
to flat spacetime intuition. This is responsible for the late 
time "tails" that are seen when a radiation packet prop- 
agates away from a black hole (cf. the late time behavior 
seen in Fig.[T]). These tails can be regarded, heuristically, 
as radiation whose propagation to large radius was de- 
layed by scattering off the spacetime. It also causes the 
integrated flux to depend on and vary with the radius at 
which it is measured. 

We now examine how our fluxes vary with respect to 
extraction radius. Tables |V] and IVII present the fluxes 
measured for four representative strong-field orbits. In 
each case, we measure E for the |m| = 2 and \m\ = 3 
modes at extraction radii R/M = 100, 200, 300, 400, 
500, and 600. These data are then fit to the ansatz 



I - q (mfloibR)' 



(5.1) 
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0.03 - 



■ta 



r 0.009 r 



0.006 
0.005 



0.003 



0.0009 r 
0.0006 



a = 0.8M, r„ = 5M 



-1 



FIG. 5: An ilustration of our code's convergence behavior. 
We show the fractional deviation in energy flux in the |m| — 2 
mode, measured ai R — 250M as a function of grid spacing. 
The grid is controlled by the integer b using 5r* = 0.0625 x 
2-*''*, S6 = Tv/30 X 2-''^*, with b G [-1,...,4]. The upper 
data set is for fluxes measured from an orbit with a = 0.8M, 
ro = 5M; the lower set is for a = 0.9M, ro = 4.64M. For each 
data set, the dotted line represents what we would expect for 
perfect second-order convergence (fit arbitrarily to the data 
for & = 0); the large dots represent our actual convergence 
data. 



The parameters g, p, and Eoa are determined by the fit. 
Notice that Eoa represents the flux that (according to 
this ansatz) would be measured infinitely far away. Note 
that this form was suggested to us by L. M. Burko [33 |. 
and replaces a previous version which used [vo-ch I RY 
rather than (mr2orb-R)~^- The two forms can be eas- 
ily related to one another; however, the form involving 
mOorb emphasizes that it is the asymptotic behavior of 
the mode, rather than a property of the orbit, that sets 
E. This form should also be more readily extendable to 
non-circular orbits. 

Figure 5 shows an example of how well this fit works 
for one of the cases given in Table fVl (rn/M = 10, 
a/M = 0.99, m = 2). Pragmatically, this ansatz ap- 
pears to fit the data quite well; the quality shown in Fig. 
iniis typical for the data that we examined. Interestingly, 
we find in all cases that the exponent p ~ 2, independent 
of black hole spin, orbit radius, or mode number. De- 
tailed calculations which we will present elsewhere [35| 
shows that this corresponds to the dominant correction 
to the asymptotic behavior of ip^. Such behavior was 
also demonstrated by Newman and Unti [3(^ (although 



our results do not currently agree with theirs when a ^ 0; 
we are investigating this discrepancy). 



X 10 ■ 



Energy flux Vs Extraction radius 



2.2 



2.15 



2.1 - 



2.05 



100 200 300 400 500 
Extraction radius (r/M) 



600 



700 



FIG. 6: A power law fit to our numerically extracted energy 
fluxes for the case ro/M = 10, a/M = 0.99, m = +2. Numer- 
ical data is indicated by the dots; the curve is the best fit we 
obtain for the ansatz given by Eq. (|5.1|l . For this case, the 
best fit parameters are q = 7.45, p = 2.06, Eoo = 2.197X lO"'^. 



We now use the fit (|5.ip to compare the extrapolated 
measured flux Eaa to frequency-domain results EpD. 
This is shown in the last column of Table IVTl In all cases, 
the error we find is less than 1%, sometimes substantially 
less. A similar fit can be performed on the angular mo- 
mentum fiux, with similar results. We conclude that the 
fit (|5.ip accounts for finite extraction radius, providing 
an accurate estimate for the fluxes that a particle radi- 
ates to infinity. Residual errors are thus much more likely 
to be true measures of numerical error in our calculation, 
and not an artifact of the extraction. 



VI. SUMMARY AND FUTURE WORK 

We have presented a simple, new technique for mod- 
eling the Dirac delta function and its derivatives on a 
finite difference grid. This technique requires that the 
source be modeled only on a handful of points on the 
grid. Our particular goal in this analysis is to model a 
pointlike source function for the time-domain Teukolsky 
equation, appropriate to describe the smaller member of 
an extreme mass ratio binary. We emphasize that our 
models for the discrete delta and its derivatives are more 
broadly applicable than just the Teukolsky equation — 
these techniques can be used in any context that requires 
modeling a sharp, delta-like function on a finite difference 
grid. 
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TABLE V: Fluxes extracted at a sequence of radii on the numerical grid. a/M is the BH spin, Ti^jM is the orbital radius and 
|m| is the azimuthal mode. Er is the flux measured at radius RM. 



vn 1 


a/M 


rtJM 


-Eioo 


-B20O 


£"300 


-E4OO 


-E5OO 


Eaoo 


2 


0.99 


4 


2.4567 X 10" 


3 


2.4681 X 10" 


3 


2.4702 X 10" 


3 


2.4709 X 10" 


3 


2.4712 X 10" 


3 


2.4714 X 10" 


3 


2 


0.99 


10 


4.1032 X 10" 


5 


4.3209 X 10" 


5 


4.3598 X 10" 


5 


4.3767 X 10" 


5 


4.3828 X 10" 


5 


4.3861 X 10" 


5 


2 


0.90 


10 


4.1729 X 10"^ 


5 


4.3930 X 10" 


5 


4.4322 X 10" 


5 


4.4456 X 10" 


5 


4.4517 X 10" 


5 


4.4550 X 10" 


5 


2 


0.00 


12 


1.9584 X 10~ 


5 


2.1256 X 10" 


5 


2.1554 X 10" 


5 


2.1654 X 10" 


5 


2.1699 X 10" 


5 


2.1723 X 10" 


5 


3 


0.99 


4 


5.8962 X 10" 


4 


5.9278 X 10" 


4 


5.9334 X 10" 


4 


5.9353 X 10" 


4 


5.9361 X 10" 


4 


5.9364 X 10" 


4 


3 


0.99 


10 


4.5778 X 10~ 


6 


4.8558 X 10" 


6 


4.9051 X 10" 


6 


4.9220 X 10" 


6 


4.9297 X 10" 


6 


4.9339 X 10" 


6 


3 


0.90 


10 


4.6791 X 10" 


6 


4.9588 X 10" 


6 


5.0085 X 10" 


6 


5.0255 X 10" 


6 


5.0333 X 10" 


6 


5.0375 X 10" 


6 


3 


0.00 


12 


1.9528 X 10" 


6 


2.1326 X 10" 


6 


2.1650 X 10" 


6 


2.1761 X 10" 


6 


2.1812 X 10" 


6 


2.1839 X 10" 


6 



TABLE VI: Best fit parameters, Boo, p, q [appearing in Eq. (|5.1|l ] for data presented in Table [V] 



|m| 


a/M 


ro/M 


Eoa 


P 


q 


Efd 


{Eoo - Efd) /Efd 


2 


0.99 


4 


2.4718 X 10"^ 


2.04 


3.40 


2.4836 X 10"^ 


-0.0048 


2 


0.99 


10 


4.3953 X lO"'' 


1.96 


2.31 


4.3948 X 10"^ 


0.0001 


2 


0.90 


10 


4.462 X lO"'^ 


2.05 


2.70 
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0.0004 



We test this approach by solving the Tcukolsky equa- 
tion for a test body in a circular, equatorial orbit of a 
Kerr black hole. Comparing with a well tested time- 
domain code that treats the orbiting body using a trun- 
cated Gaussian, we find that this new approach is ex- 
tremely fast (often by a factor of ~ 10) and accurate. 
Using a frequency-domain code as a benchmark to com- 
pare the flux of energy carried by gravitational waves, we 
find that the code which uses the discrete delta function 
is typically a factor of 2 — 5 more accurate than the Gaus- 
sian treatment most commonly used previously. This ac- 
curacy can be improved still further (at least for fluxes) 
by using a simple fit that accounts for the variation of 
the flux with the extraction radius. Combining our new 
source function with this fitting law, we find that our 
code agrees with the frequency-domain benchmark with 
errors smaller than 1% for a large fraction of parameter 
space, sometimes significantly smaller. 

Since the goal of this analysis is to contribute to the 
modeling of EMRI gravitational- wave sources, the re- 
striction to circular and equatorial orbits, though a use- 
ful, illustrative test, is not astrophysically realistic. Since 
such binaries form through scattering processes, they are 
expected to have substantial eccentricity Q , and the sec- 
ondary's orbit should have no special alignment with the 
spin axis of the large black hole. Our next analysis will 
study how well this new technique handles such orbits. 



This realistic case is substantially more difficult to treat 
than circular, equatorial orbits, since the orbiting body 
(and our discrete delta model) very rapidly crosses back 
and forth over gridpoints in both the radial and latitu- 
dinal directions. As this paper is being completed, this 
work is well underway. Early analysis points to excel- 
lent results — we get very good results even when we 
move our discrete delta model rapidly in a dynamical 
orbit. We also plan to examine wave emission from non- 
geodesic trajectories (in preparation for using the code 
to calculate wave emission from an inspiral sequence), 
and to examine plunging orbits (with a plan to examine 
waves from the transition between late inspiral to final 
plunge dzl). 

The final goal of this work will be to compute adi- 
abatic inspiral waveforms using a hybrid of frequency- 
domain and time-domain, as described in the introduc- 
tion. With a robust time-domain code for computing 
waves from nearly arbitrary physical worldlines and with 
a robust frequency-domain code capable of "mass pro- 
ducing" radiation reaction data for generic Kerr black 
hole orbits this problem should boil down to simple a 
matter of available CPU resources. Once we are in this 
state, we hope to produce waveforms efficiently enough 
that they can be used by workers looking at problems 
in LISA data analysis and waveform measurement (e.g., 
the "Mock L/S*^ Data Challenge" [H, [H, SO, IHi ) . These 
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waveforms are likely to be useful for other astrophysical 
problems, such as computing radiation recoil from both 
the slow inspiral and the dynamic plunge. Computing 
this effect in the extreme mass ratio limit may serve as 
a precision check on recent work looking at this problem 
in full numerical relativity [H, IH, HliB, El ■ 
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